library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readr)
library(readxl)
library(stringr)
library(cowplot)
##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamp
library(paletteer)
library(ggrepel)
## Warning: package 'ggrepel' was built under R version 4.3.3
all_nickel <- read.csv("../data/all_nickel_trade_data/all_nickel_2012.csv", row.names = NULL)
for (year in 2013:2024) {
new_df = read.csv(paste0("../data/all_nickel_trade_data/all_nickel_", year, ".csv"), row.names = NULL)
all_nickel <- rbind(all_nickel, new_df)
}
all_nickel <- all_nickel %>%
select(
Year = refPeriodId,
Month = refYear,
CountryCode = reporterCode,
CountryName = reporterISO,
Flow = reporterDesc,
HSCode = isOriginalClassification,
Value = fobvalue
)
all_nickel_X <- all_nickel %>%
filter(Flow == "X") %>%
group_by(Year, Month, CountryCode, CountryName) %>%
summarise(Exports = sum(Value)) %>%
ungroup()
## `summarise()` has grouped output by 'Year', 'Month', 'CountryCode'. You can
## override using the `.groups` argument.
all_nickel_yearly_actual <- read.csv("../data/all_nickel_trade_data/all_nickel_yearly.csv", row.names = NULL) %>%
select(
Year = refPeriodId,
CountryCode = reporterCode,
CountryName = reporterISO,
Flow = reporterDesc,
HSCode = isOriginalClassification,
Value = fobvalue
) %>%
filter(Flow == "X") %>%
group_by(Year, CountryCode, CountryName) %>%
summarise(Exports = sum(Value)) %>%
ungroup()
## `summarise()` has grouped output by 'Year', 'CountryCode'. You can override
## using the `.groups` argument.
nickel_price <- read.csv("../data/PNICKUSDM_from2011.csv") %>%
mutate(Date = date(observation_date)) %>%
filter(Date >= "2012-01-01") %>%
mutate(Year = as.numeric(substring(observation_date, 1, 4)))
nickel_price_yearly <- nickel_price %>%
mutate(Year = substring(observation_date, 1, 4)) %>%
group_by(Year) %>%
mutate(Year = as.numeric(Year)) %>%
summarise(MeanP = mean(PNICKUSDM))
CPI_IDN <- read_xls("../data/API_FP.CPI.TOTL_DS2_en_excel_v2_252315.xls",
sheet = "Data", range = "A4:BQ270") %>%
filter(`Country Name` == "Indonesia") %>%
select(-`Indicator Name`, -`Indicator Code`) %>%
pivot_longer(`1960`:`2024`, names_to = "Year", values_to = "CPI") %>%
mutate(Year = as.numeric(Year)) %>%
filter(Year > 2010)
USD_to_IDR_conversion <- read.csv("../data/bis_dp_search_export_20251127-202846.csv", skip = 2) %>%
mutate(Year = 2011:2024) %>%
select(Year, IDR_USD = OBS_VALUE.Value)
real_P <- CPI_IDN %>%
left_join(nickel_price, by = c("Year")) %>%
mutate(RealP = PNICKUSDM / CPI * 100) %>%
drop_na()
real_P_yearly <- CPI_IDN %>%
left_join(nickel_price_yearly, by = c("Year")) %>%
mutate(RealP = MeanP / CPI * 100) %>%
drop_na()
ggplot(real_P) +
# annotate("rect", xmin=as.Date("2014-01-01"), xmax=as.Date("2016-12-01"),
# ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
# annotate("rect", xmin=as.Date("2020-01-01"), xmax=as.Date("2024-12-31"),
# ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
geom_line(aes(x = Date, y = RealP/1e3), alpha = 0.4) +
scale_x_date(
limits = as.Date(c(NA, "2024-12-31")),
date_breaks = "1 year",
date_labels = "%Y"
) +
geom_line(
data = real_P_yearly %>%
filter(
Year >= 2012,
Year <= 2024
) %>%
mutate(Date = as.Date(paste0(Year, "-06-01"))),
aes(x = Date, y = RealP/1e3),
colour = "tomato2", size = 1.1
) +
# geom_point(
# data = nickel_price_yearly %>%
# filter(
# Year >= 2012,
# Year <= 2025
# ) %>%
# mutate(Date = as.Date(paste0(Year, "-06-01"))),
# aes(x = Date, y = MeanP/1e3),
# colour = "tomato2"
# ) +
labs(
title = "Yearly Mean Price of Nickel from 2012 to 2025",
subtitle = "Price of Nickel Futures on LME, adjusted by Indonesian CPI (100 in 2010)",
x = "Date", y = "Price (k USD per Metric Ton)",
caption = "Source: Federal Reserve Bank of St. Louis"
) +
annotate(geom = "point", x = as.Date("2022-03-01"), y = 20803.46/1e3, colour = "tomato2") +
annotate(geom = "text", x = as.Date("2022-06-01"), y = 20600/1e3, size = 3, hjust = 0,
label = "Mar 2022: \nPrice surge") +
theme_minimal() +
theme(panel.grid.minor = element_blank())
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
max(real_P$RealP)
## [1] 20803.46
# helper function
k_highest_in_period <- function(year, k) {
res = all_nickel_yearly_actual %>%
filter(Year == year) %>%
slice_max(Exports, n = k) %>%
pull(CountryCode)
return(res)
}
k_highest_Q_in_period <- function(year, k) {
res = nickel_Q %>%
filter(Year == year) %>%
slice_max(Quantity, n = k) %>%
pull(CountryCode)
return(res)
}
# plotting helpers
plot_ban_rects <- function() {
list(
annotate("rect", xmin=2013.9, xmax=2017.1, ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey"),
annotate("rect", xmin=2019.9, xmax=Inf, ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey")
)
}
plot_event_lines <- function() {
list(
geom_vline(xintercept = 2022, linetype = 2),
annotate(geom = "text", x = 2019.85, y = 7, hjust = 1, size = 2.5,
label = "2020: Indonesian export\nban starts again"),
geom_vline(xintercept = 2020, linetype = 2),
annotate(geom = "text", x = 2022.15, y = 7.85, hjust = 0, size = 2.5,
label = "2022: Price surge \nin nickel")
)
}
all_nickel_X_yearly <- all_nickel_X %>%
group_by(Year, CountryCode) %>%
summarise(Exports = sum(Exports)) %>%
ungroup()
## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.
all_nickel_X_yearly_to_plot <- all_nickel_X_yearly %>%
filter(CountryCode != "RUS") %>%
filter(CountryCode %in% k_highest_in_period(2019, 11)) %>%
left_join(
data.frame(
CountryCode = k_highest_in_period(2019, 11),
rank = seq_along(k_highest_in_period(2019, 1))
),
by = c("CountryCode")
)
all_nickel_yearly_actual_to_plot <- all_nickel_yearly_actual %>%
filter(CountryCode != "RUS") %>%
filter(CountryCode %in% k_highest_in_period(2019, 11)) %>%
left_join(
data.frame(
CountryCode = k_highest_in_period(2019, 11),
rank = seq_along(k_highest_in_period(2019, 11))
),
by = c("CountryCode")
)
k_highest_in_period(2019, 11)
## [1] "USA" "CAN" "RUS" "IDN" "DEU" "AUS" "GBR" "FIN" "NOR" "ZWE" "PHL"
# plotting
ggplot(all_nickel_X_yearly_to_plot) +
plot_ban_rects() +
plot_event_lines() +
geom_line(aes(x = Year, y = Exports/1e9, group = CountryCode, colour = rank)) +
geom_text(
data = all_nickel_X_yearly_to_plot %>%
filter(Year == 2024) %>%
select(CountryCode, LastX = Exports, rank),
aes(x = 2024.2, y = LastX/1e9, label = CountryCode),
hjust = 0, size = 2.5
) +
scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
geom_line(
data = all_nickel_X_yearly_to_plot %>% filter(CountryCode == "IDN"),
aes(x = Year, y = Exports/1e9, group = CountryCode),
colour = "tomato", linewidth = 1
) +
geom_text(
data = all_nickel_X_yearly_to_plot %>%
filter(Year == 2024) %>%
select(CountryCode, LastX = Exports, rank) %>%
filter(CountryCode == "IDN"),
aes(x = 2024.2, y = LastX/1e9, label = CountryCode),
colour = "tomato", hjust = 0, size = 2.5, fontface = "bold"
) +
scale_x_continuous(limits = c(2012, 2024.5), breaks = 2012:2024) +
labs(
title = "Nickel Exports of Top Exporting Countries from 2012 to 2024",
subtitle = "Top 8 Exporters in Dec 2019, excluding Russia",
x = "Year", y = "Value of Exports (Billion USD)",
caption = "Source: UN Comtrade Database"
) +
theme_minimal() +
theme(legend.position = "None")
ggplot(all_nickel_yearly_actual_to_plot) +
plot_ban_rects() +
plot_event_lines() +
geom_line(aes(x = Year, y = Exports/1e9, group = CountryCode, colour = rank)) +
geom_text(
data = all_nickel_yearly_actual_to_plot %>%
filter(Year == 2024) %>%
select(CountryCode, LastX = Exports, rank),
aes(x = 2024.2, y = LastX/1e9, label = CountryCode),
hjust = 0, size = 2.5
) +
scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
geom_line(
data = all_nickel_yearly_actual_to_plot %>% filter(CountryCode == "IDN"),
aes(x = Year, y = Exports/1e9, group = CountryCode),
colour = "tomato2", linewidth = 1
) +
geom_text(
data = all_nickel_yearly_actual_to_plot %>%
filter(Year == 2024) %>%
select(CountryCode, LastX = Exports, rank) %>%
filter(CountryCode == "IDN"),
aes(x = 2024.2, y = LastX/1e9, label = CountryCode),
colour = "tomato2", hjust = 0, size = 2.5, fontface = "bold"
) +
scale_x_continuous(limits = c(2013, 2024.5), breaks = 2013:2024) +
labs(
title = "Nickel Exports of Top Exporting Countries from 2012 to 2024",
subtitle = "Top 10 Exporters in Dec 2019, excluding Russia",
x = "Year", y = "Value of Exports (Billion USD)",
caption = "Source: UN Comtrade Database"
) +
theme_minimal() +
theme(legend.position = "None")
nickel_Q <- all_nickel_yearly_actual %>%
left_join(nickel_price_yearly, by = "Year") %>%
mutate(Quantity = Exports / Year)
to_plot_nickel_Q <- nickel_Q %>%
filter(CountryCode != "RUS") %>%
filter(CountryCode %in% k_highest_Q_in_period(2019, 9)) %>%
left_join(
data.frame(
CountryCode = k_highest_Q_in_period(2019, 9),
rank = seq_along(k_highest_Q_in_period(2019, 9))
),
by = c("CountryCode")
)
ggplot(to_plot_nickel_Q) +
list(
geom_vline(xintercept = 2022, linetype = 2),
annotate(geom = "text", x = 2022.15, y = 3.9, hjust = 0, size = 2.5,
label = "2022: Price surge \nin nickel")
) +
geom_line(aes(x = Year, y = Quantity/1e6, group = CountryCode, colour = rank)) +
geom_text(
data = to_plot_nickel_Q %>%
filter(Year == 2024) %>%
select(CountryCode, LastQ = Quantity, rank),
aes(x = 2024.2, y = LastQ/1e6, label = CountryCode),
hjust = 0, size = 2.5
) +
scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
geom_line(
data = to_plot_nickel_Q %>% filter(CountryCode == "IDN"),
aes(x = Year, y = Quantity/1e6, group = CountryCode),
colour = "tomato2", linewidth = 1
) +
geom_text(
data = to_plot_nickel_Q %>%
filter(Year == 2024) %>%
select(CountryCode, LastQ = Quantity, rank) %>%
filter(CountryCode == "IDN"),
aes(x = 2024.2, y = LastQ/1e6, label = CountryCode),
colour = "tomato2", hjust = 0, size = 2.5, fontface = "bold"
) +
scale_x_continuous(limits = c(2013, 2024.5), breaks = 2013:2024) +
labs(
title = "Annual Nickel Export Quantity From Top Exporting Countries from 2013 to 2024",
subtitle = "Top 8 Exporters in 2019, excluding Russia",
x = "Year", y = "Exports (Million Metric Tonnes)",
caption = "Source: UN Comtrade Database"
) +
theme_minimal() +
theme(legend.position = "None")
ggplot(nickel_price) +
annotate("rect", xmin=as.Date("2014-01-01"), xmax=as.Date("2016-12-01"),
ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
annotate("rect", xmin=as.Date("2020-01-01"), xmax=as.Date("2025-01-01"),
ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
geom_line(aes(x = Date, y = PNICKUSDM/1e3), colour = "royalblue") +
scale_x_date(
limits = as.Date(c(NA, "2025-01-01")),
date_breaks = "1 year",
date_labels = "%Y"
) +
geom_line(
data = nickel_price_yearly %>%
filter(
Year >= 2012,
Year <= 2025
) %>%
mutate(Date = as.Date(paste0(Year, "-05-15"))),
aes(x = Date, y = MeanP/1e3),
alpha = 0.5
) +
labs(
title = "Monthly Price of Nickel from 2012 to 2025 (k USD per Metric Ton)",
subtitle = "Price of Nickel Futures on LME",
x = "Date", y = "Price"
) +
annotate(geom = "text", x = as.Date("2022-06-01"), y = 33300/1e3, size = 2, hjust = 0,
label = "Mar 2022: \nPrice surge") +
theme_minimal()
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
ggplot(nickel_price) +
# annotate("rect", xmin=as.Date("2014-01-01"), xmax=as.Date("2016-12-01"),
# ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
# annotate("rect", xmin=as.Date("2020-01-01"), xmax=as.Date("2024-12-31"),
# ymin=-Inf, ymax=Inf, alpha=0.4, fill="darkgrey") +
geom_line(aes(x = Date, y = PNICKUSDM/1e3), alpha = 0.4) +
scale_x_date(
limits = as.Date(c(NA, "2024-12-31")),
date_breaks = "1 year",
date_labels = "%Y"
) +
geom_line(
data = nickel_price_yearly %>%
filter(
Year >= 2012,
Year <= 2024
) %>%
mutate(Date = as.Date(paste0(Year, "-06-01"))),
aes(x = Date, y = MeanP/1e3),
colour = "tomato2", size = 1.1
) +
# geom_point(
# data = nickel_price_yearly %>%
# filter(
# Year >= 2012,
# Year <= 2025
# ) %>%
# mutate(Date = as.Date(paste0(Year, "-06-01"))),
# aes(x = Date, y = MeanP/1e3),
# colour = "tomato2"
# ) +
labs(
title = "Yearly Mean Price of Nickel from 2012 to 2025",
subtitle = "Price of Nickel Futures on LME",
x = "Date", y = "Price (k USD per Metric Ton)",
caption = "Source: "
) +
annotate(geom = "point", x = as.Date("2022-03-01"), y = 33924.18/1e3, colour = "tomato2") +
annotate(geom = "text", x = as.Date("2022-06-01"), y = 33300/1e3, size = 2, hjust = 0,
label = "Mar 2022: \nPrice surge") +
theme_minimal() +
theme(panel.grid.minor = element_blank())
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
nickel_Q <- all_nickel_yearly_actual %>%
left_join(nickel_price_yearly, by = c("Year")) %>%
mutate(Quantity = Exports / MeanP)
nickel_Q_to_plot <- nickel_Q %>%
filter(CountryCode %in% k_highest_in_period(2019, 10)) %>%
left_join(
data.frame(
CountryCode = k_highest_Q_in_period(2019, 10),
rank = seq_along(k_highest_Q_in_period(2019, 10))
),
by = c("CountryCode")
)
ggplot(nickel_Q_to_plot, aes(x = Year, y = Quantity)) +
geom_line(aes(group = CountryCode, colour = rank)) +
geom_text(
data = nickel_Q_to_plot %>%
filter(Year == 2024) %>%
select(CountryCode, LastQ = Quantity, rank),
aes(x = 2024.2, y = LastQ, label = CountryCode),
hjust = 0, size = 2.5
) +
scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
scale_x_continuous(limits = c(2013, 2024.5), breaks = 2013:2024) +
theme_minimal() +
theme(legend.position = "None")
ggplot(nickel_Q %>%
filter(CountryCode == "IDN"), aes(x = MeanP, y = Quantity)) +
geom_text(aes(label = Year))
nickel_Q_monthly <- all_nickel_X %>%
left_join(nickel_price %>% mutate(
Year = as.numeric(substring(observation_date, 1, 4)),
Month = as.numeric(substring(observation_date, 6, 7))
), by = c("Year", "Month")) %>%
mutate(Quantity = Exports / PNICKUSDM)
nickel_Q_monthly_to_plot <- nickel_Q_monthly %>%
filter(CountryCode %in% k_highest_in_period(2019, 10)) %>%
left_join(
data.frame(
CountryCode = k_highest_in_period(2019, 10),
rank = seq_along(k_highest_in_period(2019, 10))
),
by = c("CountryCode")
)
ggplot(nickel_Q_monthly_to_plot, aes(x = Date, y = Quantity)) +
geom_line(aes(group = CountryCode, colour = rank)) +
scale_colour_gradient2(low = "orange", mid = "mediumseagreen", high = "royalblue", midpoint = 6) +
scale_x_date(
limits = as.Date(c(NA, "2025-01-01")),
date_breaks = "1 year",
date_labels = "%Y"
) +
theme_minimal() +
theme(legend.position = "None")
library(fixest)
## Warning: package 'fixest' was built under R version 4.3.3
reg_df <- all_nickel_yearly_actual %>%
left_join(nickel_price_yearly, by = "Year") %>%
mutate(
Quantity = Exports / MeanP,
ExportBan = ifelse((((Year >= 2014) & (Year <= 2016)) | (Year >= 2020) & (CountryCode == "IDN")), 1, 0),
COVID_dummy = ifelse(Year == 2020, 1, 0)
) %>%
mutate(
Int = ExportBan * MeanP
)
reg_df <- all_nickel_yearly_actual %>%
left_join(nickel_price_yearly, by = "Year") %>%
mutate(
Quantity = Exports / MeanP,
ExportBan = ifelse(((Year >= 2014) & (CountryCode == "IDN")), 1, 0),
COVID_dummy = ifelse(Year == 2020, 1, 0)
) %>%
mutate(
Int = ExportBan * MeanP
)
extra_df <- reg_df %>%
group_by(CountryCode) %>%
mutate(mean_Q = mean(Quantity)) %>%
ungroup() %>%
mutate(translated_Q = Quantity - mean_Q) %>%
group_by(Year) %>%
mutate(adj = (ifelse(CountryCode == "IDN", 1, 115)) ** 0.5) %>%
ungroup() %>%
group_by(Year, I_IDN = CountryCode == "IDN") %>%
summarise(mean_Q = sum(translated_Q) / mean(adj)) %>%
pivot_wider(names_from = I_IDN, values_from = mean_Q) %>%
mutate(diff = `TRUE` - `FALSE`)
## `summarise()` has grouped output by 'Year'. You can override using the
## `.groups` argument.
ggplot(reg_df, aes(x = Year, y = Quantity)) +
geom_line(aes(group = CountryCode, alpha = ifelse(CountryCode == "IDN", 1, 0.6))) +
geom_line(data = reg_df %>%
filter(CountryCode != "IDN") %>%
group_by(Year) %>%
summarise(Mean_Q = mean(Quantity)), aes(x = Year, y = Mean_Q * 10), colour = "tomato") +
geom_vline(xintercept = 2014, linetype = 2, 201) +
theme_minimal()
## Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.
ggplot(reg_df %>%
filter(CountryCode %in% k_highest_in_period(2019, 10)), aes(x = Year, y = Quantity)) +
geom_line(aes(group = CountryCode, alpha = ifelse(CountryCode == "IDN", 1, 0.6))) +
geom_line(data = reg_df %>%
filter(CountryCode != "IDN") %>%
filter(CountryCode %in% k_highest_in_period(2019, 10)) %>%
group_by(Year) %>%
summarise(Mean_Q = mean(Quantity)), aes(x = Year, y = Mean_Q), colour = "tomato") +
geom_vline(xintercept = 2014, linetype = 2, 201) +
theme_minimal()
## Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.
ggplot(reg_df %>% filter(CountryCode == "IDN"), aes(x = Year, y = log(Quantity))) +
plot_ban_rects() +
geom_line() +
geom_line(data = reg_df %>%
filter(CountryCode != "IDN") %>%
filter(CountryCode %in% k_highest_in_period(2019, 8)) %>%
group_by(Year) %>%
summarise(Mean_Q = mean(Quantity)), aes(x = Year, y = log(Mean_Q)), colour = "tomato") +
geom_vline(xintercept = 2014, linetype = 2, 201) +
scale_x_continuous(breaks = 2013:2024) +
theme_minimal()
## Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.
ggplot(extra_df, aes(x = Year, y = diff)) +
geom_line() +
geom_line(data = nickel_price_yearly, aes(x = Year, y = (MeanP - mean(MeanP)) * 20), colour = "royalblue", alpha = 0.7) +
geom_vline(xintercept = 2014, linetype = 2, 201) +
theme_minimal()
## Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.
base_feols <- feols(
Quantity ~ MeanP + ExportBan + Int + COVID_dummy | CountryCode,
data = reg_df %>% filter(Year > 2013))
## The variable 'ExportBan' has been removed because of collinearity (see $collin.var).
base_feols
## OLS estimation, Dep. Var.: Quantity
## Observations: 1,312
## Fixed-effects: CountryCode: 169
## Standard-errors: Clustered (CountryCode)
## Estimate Std. Error t value Pr(>|t|)
## MeanP -0.277131 0.108618 -2.55142 0.011621 *
## Int 14.220362 0.102693 138.47419 < 2.2e-16 ***
## COVID_dummy -2467.542032 1198.239091 -2.05931 0.041008 *
## ... 1 variable was removed because of collinearity (ExportBan)
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 15,020.6 Adj. R2: 0.866329
## Within R2: 0.143057
base_feols_log_1 <- feols(
log(Quantity) ~ log(MeanP) + ExportBan + ExportBan * log(MeanP) + COVID_dummy | CountryCode,
data = reg_df %>% filter(CountryCode %in% k_highest_in_period(2019, 1)))
## The variables 'ExportBan' and 'log(MeanP):ExportBan' have been removed because of collinearity (see $collin.var).
log_1_ts = base_feols_log_1$coefficients / base_feols_log_1$se
df_temp <- data.frame(k = 1, t_p = log_1_ts["log(MeanP)"], t_ban = log_1_ts["ExportBan"], t_int = log_1_ts["log(MeanP):ExportBan"])
for (i in 2:50) {
feols = feols(
log(Quantity) ~ log(MeanP) + ExportBan + ExportBan * log(MeanP) + COVID_dummy | CountryCode,
data = reg_df %>% filter(CountryCode %in% k_highest_in_period(2019, i)))
Ts = feols$coefficients / feols$se
new_row = data.frame(
k = i,
t_p = Ts["log(MeanP)"],
t_ban = Ts["ExportBan"],
t_int = Ts["log(MeanP):ExportBan"]
)
df_temp <- rbind(df_temp, new_row)
}
## The variables 'ExportBan' and 'log(MeanP):ExportBan' have been removed because of collinearity (see $collin.var).
## The variables 'ExportBan' and 'log(MeanP):ExportBan' have been removed because of collinearity (see $collin.var).
ggplot(df_temp, aes(x = k)) +
geom_hline(yintercept = -1.96, colour = "tomato", linetype = 2) +
geom_hline(yintercept = 1.96, colour = "tomato", linetype = 2) +
geom_line(aes(y = t_p), colour = "goldenrod") +
geom_line(aes(y = t_ban), colour = "seagreen") +
geom_line(aes(y = t_int), colour = "royalblue")
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_line()`).
reg_df_1 <- all_nickel_yearly_actual %>%
left_join(nickel_price_yearly, by = "Year") %>%
mutate(
Quantity = Exports / MeanP,
After = ifelse(((Year >= 2014) & (Year <= 2016)) | (Year >= 2020), 1, 0),
ExportBan = ifelse(CountryCode == "IDN", 1, 0),
COVID_dummy = ifelse(Year == 2020, 1, 0)
) %>%
mutate(
B1 = After,
B2 = After * ExportBan,
B3 = log(MeanP),
B4 = log(MeanP) * After,
B5 = log(MeanP) * ExportBan,
B6 = log(MeanP) * ExportBan * After,
Y = log(Quantity)
)
reg_df_1 <- all_nickel_yearly_actual %>%
left_join(nickel_price_yearly, by = "Year") %>%
filter(Year >= 2013) %>%
mutate(
Quantity = Exports / MeanP,
After = ifelse(Year >= 2020, 1, 0),
Indo = ifelse(CountryCode == "IDN", 1, 0),
COVID_dummy = ifelse(Year == 2020, 1, 0)
) %>%
mutate(
B1 = After,
B2 = Indo,
B3 = After * Indo,
B4 = log(MeanP),
B5 = log(MeanP) * After,
B6 = log(MeanP) * Indo,
B7 = log(MeanP) * Indo * After,
Y = log(Quantity)
)
base_feols_log <- feols(
Y ~ B1 + B3 + B4 + B6 + B7 | CountryCode,
data = reg_df_1)
base_feols_log
## OLS estimation, Dep. Var.: Y
## Observations: 1,428
## Fixed-effects: CountryCode: 172
## Standard-errors: Clustered (CountryCode)
## Estimate Std. Error t value Pr(>|t|)
## B1 0.084486 1.770513e-01 4.771850e-01 6.3384e-01
## B3 -8.228401 1.770513e-01 -4.647466e+01 < 2.2e-16 ***
## B4 -0.069190 1.857459e-01 -3.725000e-01 7.0998e-01
## B6 1.014670 1.857459e-01 5.462678e+00 1.6322e-07 ***
## B7 0.851569 2.840000e-14 2.999195e+13 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.61691 Adj. R2: 0.875782
## Within R2: 8.945e-4
base_feols_log <- feols(
Y ~ B1 + B3 + B4 + B5 + B6 + B7 | CountryCode,
data = reg_df_1 %>% filter(CountryCode != "RUS")
# %>% filter(CountryCode %in% k_highest_in_period(2019, 10))
, vcov = "hc1")
etable(base_feols_log)
## base_feols_log
## Dependent Var.: Y
##
## B1 -12.43** (4.326)
## B3 4.290 (14.07)
## B4 -0.6902* (0.3266)
## B5 1.296** (0.4484)
## B6 1.636. (0.8875)
## B7 -0.4448 (1.440)
## Fixed-Effects: -----------------
## CountryCode Yes
## _______________ _________________
## S.E. type Heteroskeda.-rob.
## Observations 1,419
## R2 0.89005
## Within R2 0.00690
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# robust to COVID control
base_feols_log <- feols(
Y ~ B1 + B3 + B4 + B5 + B6 + B7,
data = reg_df_1 %>% filter(CountryCode != "RUS")
, vcov = "hc1")
etable(base_feols_log)
## base_feols_log
## Dependent Var.: Y
##
## Constant 6.559 (8.473)
## B1 -6.678 (12.38)
## B3 -5.531 (13.67)
## B4 -0.2337 (0.8970)
## B5 0.7109 (1.281)
## B6 0.7494*** (0.0226)
## B7 0.5705 (1.373)
## _______________ __________________
## S.E. type Heteroskedas.-rob.
## Observations 1,419
## R2 0.01931
## Adj. R2 0.01515
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
base_feols_log <- feols(
Y ~ B3 + B4 | CountryCode + Year,
data = reg_df_1 %>% filter(!CountryCode %in% c("RUS", "PHL"))
, vcov = "hc1")
## The variable 'B4' has been removed because of collinearity (see $collin.var).
base_feols_log <- feols(
Y ~ B3 | CountryCode + Year,
data = reg_df_1 %>% filter(!CountryCode %in% c("RUS")))
etable(base_feols_log)
## base_feols_log
## Dependent Var.: Y
##
## B3 0.5419*** (0.1602)
## Fixed-Effects: ------------------
## CountryCode Yes
## Year Yes
## _______________ __________________
## S.E.: Clustered by: CountryCode
## Observations 1,419
## R2 0.89085
## Within R2 0.00023
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
base_feols_log <- feols(
Y ~ B3 + B7 | CountryCode + Year,
data = reg_df_1 %>% filter(!CountryCode %in% c("RUS")))
etable(base_feols_log)
## base_feols_log
## Dependent Var.: Y
##
## B3 -11.34*** (2.901)
## B7 1.207*** (0.2949)
## Fixed-Effects: -----------------
## CountryCode Yes
## Year Yes
## _______________ _________________
## S.E.: Clustered by: CountryCode
## Observations 1,419
## R2 0.89086
## Within R2 0.00032
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lm1 <- lm(Quantity ~ MeanP + After + I(MeanP * After) + COVID_dummy, data = reg_df_1 %>%
filter(CountryCode == "IDN") %>%
mutate(After = Year >= 2022))
summary(lm1)
##
## Call:
## lm(formula = Quantity ~ MeanP + After + I(MeanP * After) + COVID_dummy,
## data = reg_df_1 %>% filter(CountryCode == "IDN") %>% mutate(After = Year >=
## 2022))
##
## Residuals:
## Min 1Q Median 3Q Max
## -36422 -23966 -5748 11498 76238
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 65978.519 71617.087 0.921 0.387559
## MeanP 2.166 5.130 0.422 0.685539
## AfterTRUE 858224.099 158620.130 5.411 0.000997 ***
## I(MeanP * After) -29.432 8.293 -3.549 0.009357 **
## COVID_dummy -37226.724 44251.115 -0.841 0.428001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 41720 on 7 degrees of freedom
## Multiple R-squared: 0.9337, Adjusted R-squared: 0.8959
## F-statistic: 24.66 on 4 and 7 DF, p-value: 0.0003196
lm1 <- lm(log(Quantity) ~ log(MeanP) + After + I(log(MeanP) * After) + COVID_dummy, data = reg_df_1 %>%
filter(CountryCode == "IDN") %>%
mutate(After = Year >= 2022))
summary(lm1)
##
## Call:
## lm(formula = log(Quantity) ~ log(MeanP) + After + I(log(MeanP) *
## After) + COVID_dummy, data = reg_df_1 %>% filter(CountryCode ==
## "IDN") %>% mutate(After = Year >= 2022))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37025 -0.22855 -0.00354 0.05554 0.62909
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.7996 5.9376 1.314 0.230
## log(MeanP) 0.3787 0.6249 0.606 0.564
## AfterTRUE 21.7010 13.6080 1.595 0.155
## I(log(MeanP) * After) -2.0668 1.3793 -1.498 0.178
## COVID_dummy -0.4307 0.3991 -1.079 0.316
##
## Residual standard error: 0.3758 on 7 degrees of freedom
## Multiple R-squared: 0.8209, Adjusted R-squared: 0.7186
## F-statistic: 8.023 on 4 and 7 DF, p-value: 0.009411
base_feols_log <- feols(
log(Quantity) ~ log(MeanP) + ExportBan + ExportBan * log(MeanP) + COVID_dummy,
data = reg_df %>% filter(CountryCode %in% k_highest_in_period(2019, 10)),
vcov = "hc1")
base_feols_log
## OLS estimation, Dep. Var.: log(Quantity)
## Observations: 117
## Standard-errors: Heteroskedasticity-robust
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.764421 2.766217 5.698909 9.8986e-08 ***
## log(MeanP) -0.436004 0.289096 -1.508163 1.3433e-01
## ExportBan -18.079726 4.225375 -4.278845 3.9758e-05 ***
## COVID_dummy -0.030677 0.178009 -0.172335 8.6348e-01
## log(MeanP):ExportBan 1.888710 0.450143 4.195802 5.4666e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.648308 Adj. R2: 0.032696
base_feols_log$coefficients
## (Intercept) log(MeanP) ExportBan
## 15.76442127 -0.43600404 -18.07972613
## COVID_dummy log(MeanP):ExportBan
## -0.03067732 1.88870974
base_feols_log$se
## (Intercept) log(MeanP) ExportBan
## 2.7662174 0.2890960 4.2253751
## COVID_dummy log(MeanP):ExportBan
## 0.1780094 0.4501427
## attr(,"type")
## [1] "Heteroskedasticity-robust"
Reg formula:
\[ ln(Q) = \beta_{0} + \beta_{1}After_{t} + \beta_{2}After_{t} * ExportBan_{i} + \beta_{3}ln(P_{t}) + \beta_{4}ln(P_{t}) * After_{t} + \beta_{5}ln(Price_{t}) * ExportBan_{i} + \beta_{6}(all 3) \]
reg_df_3 <- all_nickel_yearly_actual %>%
left_join(real_P_yearly, by = "Year") %>%
filter(Year >= 2013) %>%
mutate(
Quantity = Exports / RealP,
After = ifelse(Year >= 2022, 1, 0),
Indo = ifelse(CountryCode == "IDN", 1, 0),
COVID_dummy = ifelse(Year == 2020, 1, 0)
) %>%
mutate(
B1 = After,
B2 = Indo,
B3 = After * Indo,
B4 = log(RealP),
B5 = log(RealP) * After,
B6 = log(RealP) * Indo,
B7 = log(RealP) * Indo * After,
Y = log(Quantity)
)
feols_log_3 <- feols(
Y ~ B4 + B5 + B6 + B7 | CountryCode,
data = reg_df_3 %>% filter(!CountryCode %in% c("RUS")))
etable(feols_log_3)
## feols_log_3
## Dependent Var.: Y
##
## B4 -0.6432** (0.2125)
## B5 0.0666*** (0.0191)
## B6 0.5283* (0.2125)
## B7 0.0978*** (0.0191)
## Fixed-Effects: ------------------
## CountryCode Yes
## _______________ __________________
## S.E.: Clustered by: CountryCode
## Observations 1,419
## R2 0.89042
## Within R2 0.02079
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
feols_log_event_study_0 <- feols(
Y ~ i(Year, ref = 2021) * Indo | CountryCode,
data = reg_df_3 %>% filter(!CountryCode %in% c("RUS")))
## The variable 'Indo' has been removed because of collinearity (see $collin.var).
etable(feols_log_event_study_0)
## feols_log_event..0
## Dependent Var.: Y
##
## Year = 2013 -0.3945 (0.2970)
## Year = 2014 -0.2402 (0.2902)
## Year = 2015 -0.2361 (0.2820)
## Year = 2016 0.2755 (0.2358)
## Year = 2017 0.1668 (0.2316)
## Year = 2018 0.0020 (0.1968)
## Year = 2019 0.1177 (0.2069)
## Year = 2020 -0.1629 (0.1827)
## Year = 2022 0.2480 (0.1574)
## Year = 2023 0.5366** (0.1880)
## Year = 2024 0.4216. (0.2491)
## Indo x Year = 2013 1.024*** (0.2970)
## Indo x Year = 2014 -0.0161 (0.2902)
## Indo x Year = 2015 0.0449 (0.2820)
## Indo x Year = 2016 -0.5227* (0.2358)
## Indo x Year = 2017 -0.1606 (0.2316)
## Indo x Year = 2018 0.3752. (0.1968)
## Indo x Year = 2019 0.5277* (0.2069)
## Indo x Year = 2020 -0.0237 (0.1827)
## Indo x Year = 2022 0.9876*** (0.1574)
## Indo x Year = 2023 1.056*** (0.1880)
## Indo x Year = 2024 1.600*** (0.2491)
## Fixed-Effects: ------------------
## CountryCode Yes
## __________________ __________________
## S.E.: Clustered by: CountryCode
## Observations 1,419
## R2 0.89116
## Within R2 0.02740
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
feols_log_event_study <- feols(
Y ~ i(Year, ref = 2021) * Indo | CountryCode + Year,
data = reg_df_3 %>% filter(!CountryCode %in% c("RUS")))
## The variables 'Year::2013', 'Year::2014', 'Year::2015', 'Year::2016', 'Year::2017', 'Year::2018' and 6 others have been removed because of collinearity (see $collin.var).
etable(feols_log_event_study)
## feols_log_event_..
## Dependent Var.: Y
##
## Indo x Year = 2013 1.024*** (0.2970)
## Indo x Year = 2014 -0.0161 (0.2902)
## Indo x Year = 2015 0.0449 (0.2820)
## Indo x Year = 2016 -0.5227* (0.2358)
## Indo x Year = 2017 -0.1606 (0.2316)
## Indo x Year = 2018 0.3752. (0.1968)
## Indo x Year = 2019 0.5277* (0.2069)
## Indo x Year = 2020 -0.0237 (0.1827)
## Indo x Year = 2022 0.9876*** (0.1574)
## Indo x Year = 2023 1.056*** (0.1880)
## Indo x Year = 2024 1.600*** (0.2491)
## Fixed-Effects: ------------------
## CountryCode Yes
## Year Yes
## __________________ __________________
## S.E.: Clustered by: CountryCode
## Observations 1,419
## R2 0.89116
## Within R2 0.00119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
event_study_plotting_data <- data.frame(
Coefficient = feols_log_event_study$coefficients,
se = feols_log_event_study_0$se[1:11]
) %>%
mutate(Year = c(2013:2020, 2022:2024)) %>%
mutate(p95_length = se * 1.96)
ref_row <- data.frame(
Year = 2021,
Coefficient = 0,
se = 0,
p95_length = 0
)
plot_df <- rbind(event_study_plotting_data, ref_row)
ggplot(plot_df, aes(
x = Year,
y = Coefficient,
ymin = Coefficient - p95_length,
ymax = Coefficient + p95_length
)) +
geom_hline(yintercept = 0, colour = "grey") +
geom_vline(xintercept = 2021, linetype = 2) +
geom_pointrange() +
geom_line() +
scale_x_continuous(breaks = 2013:2024) +
labs(title = "Event Study of Effect of Downstreaming Policies on Export Level") +
theme_minimal()
feols_log_event_study_2 <- feols(
Y ~ i(Year, ref = 2021) * Indo * log(MeanP) | CountryCode,
data = reg_df_1 %>% filter(!CountryCode %in% c("RUS")))
## The variables 'Indo', 'log(MeanP)', 'log(MeanP):Year::2013', 'log(MeanP):Year::2015', 'log(MeanP):Year::2016', 'log(MeanP):Year::2017' and 18 others have been removed because of collinearity (see $collin.var).
etable(feols_log_event_study_2)
## feols_log_event..2
## Dependent Var.: Y
##
## Year = 2013 -0.1029 (0.2972)
## Year = 2014 -0.0104 (0.2518)
## Year = 2015 -0.0683 (0.2821)
## Year = 2016 0.4087. (0.2359)
## Year = 2017 0.2627 (0.2316)
## Year = 2018 0.0664 (0.1968)
## Year = 2019 0.1522 (0.2070)
## Year = 2020 -0.1474 (0.1827)
## Year = 2022 0.2068 (0.1574)
## Year = 2023 0.4593* (0.1881)
## Year = 2024 0.3228 (0.2492)
## Indo x Year = 2013 1.024*** (0.2972)
## Indo x Year = 2014 -0.0161 (0.2903)
## Indo x Year = 2015 0.0449 (0.2821)
## Indo x Year = 2016 -0.5227* (0.2359)
## Indo x Year = 2017 -0.1606 (0.2316)
## Indo x Year = 2018 0.3752. (0.1968)
## Indo x Year = 2019 0.5277* (0.2070)
## Indo x Year = 2020 -0.0237 (0.1827)
## Indo x Year = 2022 0.9876*** (0.1574)
## Indo x Year = 2023 1.056*** (0.1881)
## Indo x Year = 2024 1.600*** (0.2492)
## log(MeanP) x Year = 2014 -3.24e-5 (0.0045)
## Fixed-Effects: ------------------
## CountryCode Yes
## ________________________ __________________
## S.E.: Clustered by: CountryCode
## Observations 1,419
## R2 0.89096
## Within R2 0.01509
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
event_study_plotting_data <- data.frame(
Coefficient = feols_log_event_study$coefficients[12:22],
se = feols_log_event_study$se[1:11]
) %>%
mutate(Year = c(2013:2020, 2022:2024)) %>%
mutate(p95_length = se * 1.96)
ref_row <- data.frame(
Year = 2021,
Coefficient = 0,
se = 0,
p95_length = 0
)
plot_df <- rbind(event_study_plotting_data, ref_row)
ggplot(plot_df, aes(
x = Year,
y = Coefficient,
ymin = Coefficient - p95_length,
ymax = Coefficient + p95_length
)) +
geom_hline(yintercept = 0, colour = "grey") +
geom_vline(xintercept = 2021, linetype = 2) +
geom_pointrange() +
geom_line() +
scale_x_continuous(breaks = 2013:2024) +
theme_minimal()
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_pointrange()`).
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_line()`).
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
reg_df_2 <- all_nickel_yearly_actual %>%
left_join(real_P_yearly, by = "Year") %>%
filter(Year >= 2013) %>%
group_by(CountryCode) %>%
mutate(
Qt = Exports / RealP,
Pt = RealP,
) %>%
mutate(
Qt_1 = lag(Qt),
Pt_1 = lag(Pt)
) %>%
ungroup() %>%
mutate(
Indo = ifelse(CountryCode == "IDN", 1, 0),
Downstream = ifelse(Year >= 2020, 1, 0),
After = ifelse(Year >= 2022, 1, 0),
COVID_dummy = ifelse(Year == 2020, 1, 0),
Y = log(Qt / Qt_1),
dP = log(Pt / Pt_1)
) %>%
drop_na()
feols_log_2 <- feols(
Y ~ Indo * Downstream * dP | CountryCode,
data = reg_df_2 %>% filter(!CountryCode %in% c("RUS")) %>% filter(CountryCode %in% k_highest_in_period(2019, 20)))
## The variable 'Indo' has been removed because of collinearity (see $collin.var).
etable(feols_log_2)
## feols_log_2
## Dependent Var.: Y
##
## Downstream -0.1341. (0.0696)
## dP 0.6848 (0.7234)
## Indo x Downstream 0.3907*** (0.0696)
## Indo x dP -0.5617 (0.7234)
## Downstream x dP -0.9541 (0.7666)
## Indo x Downstream x dP 1.644* (0.7666)
## Fixed-Effects: ------------------
## CountryCode Yes
## ______________________ __________________
## S.E.: Clustered by: CountryCode
## Observations 201
## R2 0.00718
## Within R2 0.00591
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
wald(feols_log_2, "dP + dP:Indo = 0")
## [1] NA
https://pubs.usgs.gov/myb/vol3/2023/myb3-2023-indonesia.pdf
num_smelters <- data.frame(
Year = c(2013, 2014, 2020, 2023),
completed = c(2, 3, 13, 44),
constructing = c(0, 0, 0, 44+28),
planned = c(0, 0, 0, 44+28+24)
)
num_complete <- num_smelters %>% select(Year, completed)
for (i in 2015:2024) {
if (i %in% num_smelters$Year) {next}
else if (i < 2020) {
num_complete <- rbind(
num_complete,
data.frame(
Year = i,
completed = num_smelters[2, "completed"] + (num_smelters[3, "completed"] - num_smelters[2, "completed"]) * (i - 2014) / 6
)
)
}
else if (i < 2023) {
num_complete <- rbind(
num_complete,
data.frame(
Year = i,
completed = num_smelters[3, "completed"] + (num_smelters[4, "completed"] - num_smelters[3, "completed"]) * (i - 2020) / 3
)
)
}
else {
num_complete <- rbind(
num_complete,
data.frame(
Year = i,
completed = num_smelters[4, "completed"] + (num_smelters[4, "completed"] - num_smelters[3, "completed"]) * (i - 2023) / 3
)
)
}
}
ggplot(num_smelters, aes(x = Year)) +
geom_col(aes(y = completed), fill = "royalblue") +
geom_point(aes(y = completed)) +
geom_line(aes(y = completed)) +
geom_text(aes(y = completed + 3, label = completed)) +
labs(title = "Number of Nickel Smelters in Indonesia",
x = "Year", y = "Number of Smelters",
caption = "Source: Chung, 2023") +
scale_x_continuous(breaks = 2013:2024) +
theme_minimal()
ggplot(num_complete, aes(x = Year)) +
geom_col(aes(y = completed), fill = "royalblue") +
geom_point(aes(y = completed)) +
geom_line(aes(y = completed)) +
geom_text(aes(y = completed + 3, label = completed)) +
labs(title = "Number of Nickel Smelters in Indonesia",
x = "Year", y = "Number of Smelters",
caption = "Source: ") +
scale_x_continuous(breaks = 2013:2024) +
theme_minimal()
reg_df_4 <- all_nickel_yearly_actual %>%
left_join(real_P_yearly, by = "Year") %>%
filter(Year >= 2014) %>%
mutate(
Quantity = Exports / RealP,
Downstream = ifelse(Year >= 2020, 1, 0),
Surge = ifelse(Year %in% c(2022), 1, 0),
Decline = ifelse(Year %in% c(2023, 2024), 1, 0),
After = ifelse(Year >= 2022, 1, 0),
Indo = ifelse(CountryCode == "IDN", 1, 0),
COVID_dummy = ifelse(Year == 2020, 1, 0)
) %>%
left_join(num_complete, by = "Year")
feols_log_4 <- feols(
log(Quantity) ~ Surge * log(RealP) + Decline * log(RealP) + Indo * Downstream | CountryCode,
data = reg_df_4
)
## The variables 'Indo' and 'Surge:log(RealP)' have been removed because of collinearity (see $collin.var).
etable(feols_log_4)
## feols_log_4
## Dependent Var.: log(Quantity)
##
## Surge 0.5863** (0.1785)
## log(RealP) -0.5601* (0.2550)
## Decline -8.067 (7.006)
## Downstream -0.0347 (0.1664)
## log(RealP) x Decline 0.9324 (0.7459)
## Indo x Downstream 0.6899*** (0.1503)
## Fixed-Effects: ------------------
## CountryCode Yes
## ____________________ __________________
## S.E.: Clustered by: CountryCode
## Observations 1,312
## R2 0.89721
## Within R2 0.01864
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
feols_log_4 <- feols(
log(Quantity) ~ Indo * After * log(RealP) + Indo * Downstream | CountryCode,
data = reg_df_4
)
## The variable 'Indo' has been removed because of collinearity (see $collin.var).
etable(feols_log_4)
## feols_log_4
## Dependent Var.: log(Quantity)
##
## After -1.408 (5.411)
## log(RealP) -0.5644* (0.2572)
## Downstream -0.0277 (0.1668)
## Indo x After 18.16*** (5.411)
## Indo x log(RealP) 0.5496* (0.2572)
## After x log(RealP) 0.2130 (0.5705)
## Indo x Downstream -0.1192 (0.1668)
## Indo x After x log(RealP) -1.807** (0.5705)
## Fixed-Effects: -----------------
## CountryCode Yes
## _________________________ _________________
## S.E.: Clustered by: CountryCode
## Observations 1,312
## R2 0.89718
## Within R2 0.01834
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
feols_log_4 <- feols(
log(Quantity) ~ After * log(RealP) + completed,
data = reg_df_4 %>% filter(CountryCode == "IDN")
)
etable(feols_log_4)
## feols_log_4
## Dependent Var.: log(Quantity)
##
## Constant 13.27* (5.377)
## After 8.402 (14.64)
## log(RealP) -0.1980 (0.5929)
## completed 0.0152 (0.0206)
## After x log(RealP) -0.7696 (1.511)
## __________________ ________________
## S.E. type IID
## Observations 11
## R2 0.89631
## Adj. R2 0.82718
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
feols_log_4 <- feols(
log(Quantity) ~ After * log(RealP) + completed * Indo + COVID_dummy | CountryCode,
data = reg_df_4
)
## The variable 'Indo' has been removed because of collinearity (see $collin.var).
etable(feols_log_4)
## feols_log_4
## Dependent Var.: log(Quantity)
##
## After -7.554 (9.188)
## log(RealP) -0.6719* (0.2638)
## completed 0.0117 (0.0140)
## COVID_dummy -0.2516. (0.1410)
## After x log(RealP) 0.8202 (0.9386)
## completed x Indo 0.0321*** (0.0048)
## Fixed-Effects: ------------------
## CountryCode Yes
## __________________ __________________
## S.E.: Clustered by: CountryCode
## Observations 1,312
## R2 0.89744
## Within R2 0.02087
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
feols_log_4 <- feols(
log(Quantity) ~ log(RealP) + Indo * Downstream | CountryCode,
data = reg_df_4
)
## The variable 'Indo' has been removed because of collinearity (see $collin.var).
etable(feols_log_4)
## feols_log_4
## Dependent Var.: log(Quantity)
##
## log(RealP) -0.2030 (0.1972)
## Downstream 0.2311 (0.1657)
## Indo x Downstream 0.6981*** (0.1493)
## Fixed-Effects: ------------------
## CountryCode Yes
## _________________ __________________
## S.E.: Clustered by: CountryCode
## Observations 1,312
## R2 0.89570
## Within R2 0.00425
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# event study
feols_log_4_event_study <- feols(
log(Quantity) ~ log(RealP) * i(Year, ref = 2021) + Indo * i(Year, ref = 2021) | CountryCode,
data = reg_df_4
)
## The variables 'Year::2024', 'Indo', 'log(RealP):Year::2016', 'log(RealP):Year::2018', 'log(RealP):Year::2020', 'log(RealP):Year::2022' and 2 others have been removed because of collinearity (see $collin.var).
etable(feols_log_4_event_study)
## feols_log_4_even..
## Dependent Var.: log(Quantity)
##
## log(RealP) -2.230. (1.279)
## Year = 2014 0.2318 (0.2989)
## Year = 2015 -0.5113* (0.2210)
## Year = 2016 -0.8807 (0.6358)
## Year = 2017 -0.4850 (0.3922)
## Year = 2018 -0.6164. (0.3578)
## Year = 2019 -0.2812 (0.2000)
## Year = 2020 -0.7795* (0.3288)
## Year = 2022 0.9138. (0.4637)
## Year = 2023 0.7071** (0.2608)
## log(RealP) x Year = 2014 -0.0158 (0.0121)
## log(RealP) x Year = 2015 -0.0351* (0.0150)
## log(RealP) x Year = 2017 -0.0452 (0.0276)
## log(RealP) x Year = 2019 -0.0161 (0.0116)
## Indo x Year = 2014 -0.0242 (0.2891)
## Indo x Year = 2015 0.0288 (0.2804)
## Indo x Year = 2016 -0.5293* (0.2318)
## Indo x Year = 2017 -0.1714 (0.2267)
## Indo x Year = 2018 0.3739. (0.1904)
## Indo x Year = 2019 0.5193* (0.2029)
## Indo x Year = 2020 -0.0236 (0.1814)
## Indo x Year = 2022 0.9813*** (0.1562)
## Indo x Year = 2023 1.055*** (0.1848)
## Indo x Year = 2024 1.592*** (0.2463)
## Fixed-Effects: ------------------
## CountryCode Yes
## ________________________ __________________
## S.E.: Clustered by: CountryCode
## Observations 1,312
## R2 0.89782
## Within R2 0.02444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
event_study_plotting_data <- data.frame(
Coefficient = feols_log_4_event_study$coefficients[12:22],
se = feols_log_4_event_study$se[1:11]
) %>%
mutate(Year = c(2013:2020, 2022:2024)) %>%
mutate(p95_length = se * 1.96)
ref_row <- data.frame(
Year = 2021,
Coefficient = 0,
se = 0,
p95_length = 0
)
plot_df <- rbind(event_study_plotting_data, ref_row)
ggplot(plot_df, aes(
x = Year,
y = Coefficient,
ymin = Coefficient - p95_length,
ymax = Coefficient + p95_length
)) +
geom_hline(yintercept = 0, colour = "grey") +
geom_vline(xintercept = 2021, linetype = 2) +
geom_pointrange() +
geom_line() +
scale_x_continuous(breaks = 2013:2024) +
theme_minimal()
1,532 trillion
pop_mult = 1.008
df_temp_pop <- data.frame(
year = c(2019, 2020, 2021, 2022, 2023, 2024, 2025, 2026, 2027, 2028, 2029, 2030, 2031, 2032, 2033, 2034, 2035),
resident = c(4026209, 4044210, 3986842, 4073239, 4149253, 4180868, 4204515, 4204515 * pop_mult,
4204515 * pop_mult ** 2, 4204515 * pop_mult ** 3, 4204515 * pop_mult ** 4, 4204515 * pop_mult ** 5, 4204515 * pop_mult ** 6, 4204515 * pop_mult ** 7, 4204515 * pop_mult ** 8, 4204515 * pop_mult ** 9, 4204515 * pop_mult ** 10),
population = c(5703569, 5685807, 5453566, 5637022, 5917648, 6036860, 6111175, 6111175 * pop_mult,
6111175 * pop_mult ** 2, 6111175 * pop_mult ** 3, 6111175 * pop_mult ** 4, 6111175 * pop_mult ** 5, 6111175 * pop_mult ** 6, 6111175 * pop_mult ** 7, 6111175 * pop_mult ** 8, 6111175 * pop_mult ** 9, 6111175 * pop_mult ** 10)
)
ggplot(df_temp_pop, aes(x = year, y = population)) +
geom_line(alpha = 0.5, colour = "royalblue", size = 1.2) +
geom_point(alpha = 0.5, colour = "royalblue", size = 3) +
geom_line(data = df_temp_pop %>% filter(year <= 2025), colour = "royalblue", size = 1.2) +
geom_point(data = df_temp_pop %>% filter(year <= 2025), colour = "royalblue", size = 3) +
theme_minimal()
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.